classdef v1_population <population
    properties
        N = 100;
        v1_neurons
    end
    properties(Transient = true)
        locations
        radii
        handles
    end
    methods
        function p = v1_population(N,neurons,type)
            if ~exist('type','var')||isempty(type)
                type = 'v1_population';
            end
            p = p@population(type);
            
            if (mod(N,1)==0)
                p.N = N;
            else
                error('v1_population:wrongType','N has to be a whole number only');
            end
            
            if isvector(neurons) && isa(neurons,'v1_neuron')
                p.v1_neurons = neurons;
            elseif isstruct(neurons) && ismember(neurons.method,{'AllRandom'})
                fn = str2func(neurons.method);
                p = fn(p,neurons.params);
            end
        end
        function p = AllRandom(p,params)
            neurons(1,p.N) = v1_neuron;
            for i = 1:p.N
                bgRate = params.bgRate(1)+rand*(params.bgRate(2)-params.bgRate(1));
                maxRate = params.maxRate(1)+rand*(params.maxRate(2)-params.maxRate(1));
                
                sfDetails.calcMethod = params.sfDetails.calcMethod;
                sfDetails.tuningParams.optSF = params.sfDetails.tuningParams.optSF(1)+ rand*(params.sfDetails.tuningParams.optSF(2)-params.sfDetails.tuningParams.optSF(1));
                sfDetails.tuningParams.sigmaSF = params.sfDetails.tuningParams.sigmaSF(1)+ rand*(params.sfDetails.tuningParams.sigmaSF(2)-params.sfDetails.tuningParams.sigmaSF(1));
                sfObj = sfTuning(sfDetails.calcMethod,sfDetails.tuningParams);
                
                tfDetails.calcMethod = params.tfDetails.calcMethod;
                tfDetails.tuningParams.optTF = params.tfDetails.tuningParams.optTF(1)+ rand*(params.tfDetails.tuningParams.optTF(2)-params.tfDetails.tuningParams.optTF(1));
                tfDetails.tuningParams.sigmaTF = params.tfDetails.tuningParams.sigmaTF(1)+ rand*(params.tfDetails.tuningParams.sigmaTF(2)-params.tfDetails.tuningParams.sigmaTF(1));
                tfObj = tfTuning(tfDetails.calcMethod,tfDetails.tuningParams);
                
                orDetails.calcMethod = params.orDetails.calcMethod;
                orDetails.tuningParams.alpha = 0:0.01:2*pi;
                orDetails.tuningParams.thetahat = params.orDetails.tuningParams.thetahat(1)+rand*(params.orDetails.tuningParams.thetahat(2)-params.orDetails.tuningParams.thetahat(1));
                orDetails.tuningParams.kappa = params.orDetails.tuningParams.kappa(1)+rand*(params.orDetails.tuningParams.kappa(2)-params.orDetails.tuningParams.kappa(1));
                orDetails.tuningParams.DSI = params.orDetails.tuningParams.DSI(1)+rand*(params.orDetails.tuningParams.DSI(2)-params.orDetails.tuningParams.DSI(1));
                orObj = orTuning(orDetails.calcMethod,orDetails.tuningParams);
                
                ctrDetails.calcMethod = params.ctrDetails.calcMethod;
                ctrDetails.tuningParams.n = params.ctrDetails.tuningParams.n(1)+rand*(params.ctrDetails.tuningParams.n(2)-params.ctrDetails.tuningParams.n(1));
                ctrDetails.tuningParams.c50 = params.ctrDetails.tuningParams.c50(1)+rand*(params.ctrDetails.tuningParams.c50(2)-params.ctrDetails.tuningParams.c50(1));
                ctrObj = ctrTuning(ctrDetails.calcMethod,ctrDetails.tuningParams);
                
                rfCen = [NaN NaN];
                rfSize = [NaN NaN];
                
                firingStats = 'poissrnd';
                
                neurons(i) = v1_neuron(bgRate,maxRate,sfObj,tfObj,ctrObj,orObj,rfCen,rfSize,firingStats);
            end
            p.v1_neurons = neurons;
        end
        function p = plot(p,params)
            if ~isfield(p.handles,'figure')||isempty(p.handles.figure)||~ishandle(p.handles.figure) % either field doesnt exist, is empty, or has been closed
                p.handles.figure = figure;
            else
                figure(p.handles.figure);
            end
            
            if ~isfield(p.handles,'axes')||isempty(p.handles.axes)||~ishandle(p.handles.axes) % either field doesnt exist, is empty, or has been closed
                p.handles.axes = axes;
            else
                axes(p.handles.axes);
            end
            
            if ~isfield(p.handles,'neurons')||isempty(p.handles.neurons)||~any(ishandle(p.handles.neurons)) % either field doesnt exist, is empty, or has been closed
                % arrange neurons and then initiate the axes;
                in.N = p.N;
                in.rNeuron = 3;
                in.fov.x = [-100 100];
                in.fov.y = [-100 100];
                out = arrangeNeurons(in,false);
                p.radii = repmat(3,1,p.N);
                p.locations = out.cenXY;
                p.handles.neurons = ellipse(p.radii,p.radii,zeros(p.N,1),p.locations(:,1),p.locations(:,2),'k');
                axis tight
            end
            
            if ~exist('params','var')||isempty(params)
                params.stim = grayScreen;
            end
            
            if ~isfield(params,'activity')
                params.activity = sample(p,params.stim);
            elseif length(params.activity)~=p.N
                error('v1_population:wrongLength','numNeurons:%d; numActivity provided:%d',p.N,length(params.activity));
            end
            
            normalizedActivity = params.activity./getMaxRates(p);
            alpha = 0:0.01:2*pi;
            beachBallColor = jet(length(alpha));
            ORs = getOptORs(p)';
            colScale = nan(p.N,4);
            
            [~, index] = min(abs(repmat(alpha,length(ORs),1)-repmat(ORs,1,length(alpha))),[],2);
            colScale(:,1:3) = beachBallColor(index,:);
            colScale(:,4) = normalizedActivity;
            
            XData = get(p.handles.neurons,'XData');
            YData = get(p.handles.neurons,'YData');
            axes(p.handles.axes); hold on;
            
            if ~isfield(p.handles,'fills')||isempty(p.handles.fills)||~any(ishandle(p.handles.fills)) % either field doesnt exist, is empty, or has been closed
                p.handles.fills(1,p.N) = nan;
                for i = 1:p.N
                    if colScale(i,4)>1
                        colScale(i,4) = 1;
                    end
                    p.handles.fills(i) = fill(XData{i},YData{i},colScale(i,1:3));set(p.handles.fills(i),'FaceAlpha',colScale(i,4));
                end
            else
                for i = 1:p.N
                    if colScale(i,4)>1
                        colScale(i,4) = 1;
                    end
                    set(p.handles.fills(i),'FaceAlpha',colScale(i,4));
                end
            end
            
        end
        function out = sample(p,stimulus)
            out = nan(1,p.N);
            for i = 1:p.N
                out(i) = sample(p.v1_neurons(i),stimulus);
            end
        end
        function out = getMaxRates(p)
            out = nan(1,p.N);
            for i = 1:p.N
                out(i) = p.v1_neurons(i).maxRate;
            end
        end
        function out = getOptORs(p)
            out = nan(1,p.N);
            for i = 1:p.N
                out(i) = p.v1_neurons(i).orTuning.getOptOr;
            end
        end
    end
end